function NotchBeam1
clc
%close all
addpath('../');
addpath('../../');
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% create mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mesh = fem2.Mesh();
n1	=	mesh.addNode	(	fem2.Node	(	0	,	0	)	)	;
n2	=	mesh.addNode	(	fem2.Node	(	0.95	,	0	)	)	;
n3	=	mesh.addNode	(	fem2.Node	(	0.95	,	0.15	)	)	;
n4	=	mesh.addNode	(	fem2.Node	(	1.05	,	0.15	)	)	;
n5	=	mesh.addNode	(	fem2.Node	(	1.05	,	0	)	)	;
n6	=	mesh.addNode	(	fem2.Node	(	2	,	0	)	)	;
n7	=	mesh.addNode	(	fem2.Node	(	2	,	1	)	)	;
n8	=	mesh.addNode	(	fem2.Node	(	1.05	,	1	)	)	;
n9	=	mesh.addNode	(	fem2.Node	(	1.05	,	0.85	)	)	;
n10	=	mesh.addNode	(	fem2.Node	(	0.95	,	0.85	)	)	;
n11	=	mesh.addNode	(	fem2.Node	(	0.95	,	1	)	)	;
n12	=	mesh.addNode	(	fem2.Node	(	0	,	1	)	)	;
n14	=	mesh.addNode	(	fem2.Node	(	0	,	0.15	)	)	;
n15	=	mesh.addNode	(	fem2.Node	(	2	,	0.85	)	)	;
n18	=	mesh.addNode	(	fem2.Node	(	0	,	0.85	)	)	;
n20	=	mesh.addNode	(	fem2.Node	(	2	,	0.15	)	)	;
mesh.addFace(fem2.Face([n1 n2 n3]));
mesh.addFace(fem2.Face([n1 n3 n14]));
mesh.addFace(fem2.Face([n14 n3 n10]));
mesh.addFace(fem2.Face([n14 n10 n18]));
mesh.addFace(fem2.Face([n18 n10 n11]));
mesh.addFace(fem2.Face([n18 n11 n12]));
mesh.addFace(fem2.Face([n3 n4 n9]));
mesh.addFace(fem2.Face([n3 n9 n10]));
mesh.addFace(fem2.Face([n5 n6 n20]));
mesh.addFace(fem2.Face([n5 n20 n4]));
mesh.addFace(fem2.Face([n4 n20 n15]));
mesh.addFace(fem2.Face([n4 n15 n9]));
mesh.addFace(fem2.Face([n9 n15 n7]));
mesh.addFace(fem2.Face([n9 n7 n8]));

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% config for mesh
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mesh.setNumberOfFields(2);
%set id for nodes
for i=1:mesh.countNodes()
    mesh.getNode(i-1).setId(i);
end
%set id for faces
for i=1:mesh.countFaces()
    mesh.getFace(i-1).setId(i);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% conditions
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
du = 0.65;
lambda = 0:0.1:1.1;
meshUtil = fem2.MeshUtilities(mesh);
c1 = fem2.Constraint([0 0]);
c2 = fem2.Constraint([0 1]);
c3 = fem2.Constraint([0 1]);
c3.setValue(0,du);
edge1 = meshUtil.seekNodesOnSurface(1,0,0);
edge2 = meshUtil.seekNodesOnSurface(1,0,-2);
mesh.setConstraint(c2,edge1);
meshUtil.seekNode([0 0]).setConstraint(c1);
mesh.setConstraint(c3,edge2);
% plotMesh(mesh);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% element declaration
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
M = fem2.Model(mesh);
nf = mesh.countFaces();
E = 2;
nu = 0.3;
H = -1;
ft = 0.91;
C = C_plane_stress(E,nu);
for i=1:nf
    e = fem2.MatlabStructuralElement(mesh.getFace(i-1));
    M.addElement(e);
    %create struct to contain element state variables
    es = createElementContainer(C,H,ft);
    elements(i) = es;
end
%deactivate and activate elements
% for i=1:nf
%     c = computeFaceCenter(mesh.getFace(i-1));
%     if (c(1)>=0 && c(1)<=0.6 && c(2)>=-0.6 && c(2)<=-0.4)
%         elements(i).isActivated = 1;
%     else
%         elements(i).isActivated = 0;
%     end
% end

for i=1:edge1.size()
    edge1nodes(i) = edge1.get(i-1).getId();
end

postProcessInfo.infoList = {'resistant force'};
postProcessInfo.infoParam = {edge1nodes};

[elements,output] = solving_strategy_static_2(elements,M,...
    lambda,50,1e-6,postProcessInfo);

figure
res_node_1(:,:) = -output.resistant_force(:,1,:)
res_node_3(:,:) = -output.resistant_force(:,3,:)
delta = lambda*du;
plot(delta,res_node_3);

plotMeshWithCrack(M,elements,{'trimesh','conditions','id'});

rmpath('../');
rmpath('../../');

